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ABSTRACT 

We use weak gravitational lensing to measure mean mass profiles around Locally Brightest 
Galaxies (LBGs). These are selected from the SDSS/DR7 spectroscopic and photometric cata¬ 
logues to be brighter than any neighbour projected within 1.0 Mpc and differing in redshift by 
< 1000 km/s. Most (> 83%) are expected to be the central galaxies of their dark matter halos. 
Previous stacking analyses have used this LBG sample to measure mean Sunyaev-Zeldovich 
flux and mean X-ray luminosity as a function of LBG stellar mass. In both cases, a simulation 
of the formation of the galaxy population was used to estimate effective halo mass for LBGs 
of given stellar mass, allowing the derivation of scaling relations between the gas properties 
of halos and their mass. By comparing results from a variety of simulations to our lensing 
data, we show that this procedure has significant model dependence reflecting: (i) the failure 
of any given simulation to reproduce observed galaxy abundances exactly; (ii) a dependence 
on the cosmology underlying the simulation; and (iii) a dependence on the details of how 
galaxies populate halos. We use our lensing results to recalibrate the scaling relations, elimi¬ 
nating most of this model dependence and explicitly accounting both for residual modelling 
uncertainties and for observational uncertainties in the lensing results. The resulting scaling 
relations link the mean gas properties of dark halos to their mass over an unprecedentedly 
wide range, 10^^ ® < Msoo/Mq < 10^^'^, and should fairly and robustly represent the full 
halo population. 


1 INTRODUCTION 

Recent measurements of fluctuations in the microwave back¬ 
ground have determined the fraction of all matter that is bary- 
on ic to high accuracy, flt/fl m = 0.158 ± 0.002, according 
to IPlanck CollaborationI ( l2015h . At low redshift, the baryon frac¬ 
tion of the richest galaxy clusters is close to this value, with 
~85% of the material in the form o f X-ray emitting gas and 
only ~15% in the form of stars (e.g. I Allen et ^l2008h . In less 
massive systems, the detected baryon fraction is lower and is 
dominated, for individual galaxy halos, by the stars in the cen¬ 
tral galaxy. This contribution maximises at about 4% within the 
virial radius for galaxies similar in mass to the Milky Way ( e.g. 
iMandelbaum et al1l2006l : Iwang & Jing|l20ld : IGuo et S]|2010l) . In 

the low redshift universe as a whole, stars account for only a 
few percent of the expected baryons jLi & Whitel l2009l) . Most 
remain undetected and are thought to be in the intergalactic 
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medium, either photoionised by the metagalactic UV background, 
or heated, enriched and ejected from galaxy halos by feedback from 
star formation and / or AGN activity (e . g. ICen & Ostrikeil [19^: 
Menard_et^ 20K; Menard & Fukugital2012l: Vogelsberger et al.l 


20l4ISchave et al.ll2015h.' 

Recent stacking resul ts based on the Planck sky maps 
I PlmckCofl^OTati^ 1201 3h and on the Rosat All Sky Survey 
I Anderson et al. 20151) have substantially extended the halo mass 
range over which signals from associated hot gas have been de¬ 
tected. Both studies stacked signals as a function of stellar mass 
around a sample of Locally Brightest Galaxies (LBGs) taken 
from the Se venth Data Release o f the Sloan Digital Sky Survey 
(SDSS/DR7 lAbazaiian et alll2009l) . An independent stacking anal¬ 
ysis of Planck data by Greco et ^ ( 1201 Sh recently confirmed the 
original estimates of mean SZ signal as a function of LBG stellar 
mass. These galaxies are selected to be brighter than any of their 
neighbours and are predominantly the central galaxies of their dark 
halos. 

The measured Sunyaev-Zeldovich (SZ) and X-ray signals as 
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a function of LBG stellar mass can be converted into scaling re¬ 
lations between halo mass and halo gas properties provided that 
effec tive halo mass is known as a fu nction of LBG stellar mass. 
Both IPlanck CollaboratiorJ (1201 3h and I Anderson et all (1201 Sh esti¬ 
mated this relation by forward modelling of sample selection and 
signal m easurement using a simulation of galaxy population evolu¬ 
tion from iGuo et alJi2013h . This simulation was tuned to reproduce 
the observed stellar mass function in the SDSS, so this procedure 
amounts to an abundance-matching calibration of halo mass. Both 
scaling relations were found to be unbroken power-laws over the 
range 10 ^^ ® < Mha(o/M 0 < but whereas the relation 

for SZ flux was consistent with detection of the cosmic baryon 
fraction in hot gas at all halo masses (e.g. IPlanck Collaborationl 
l201ll : lMarrone et al.l2012h . that for the X-ray flux was su bstantially 
steeper than self-similar behaviour would require (e.g. Dai_e^l.l 
|200^: iRvkoff et alj|2008l: IVikhlinin et Zll2009t IWang et alj|20l4 . 
suggesting that the concentration of baryons within halos must 
decrease with decreasing halo mass, i.e., the hot gas distribution 
within lower mass halos is less centrally concentrated than those 
in massive clusters. iLe Brun et al.l j2015h showed that when reso¬ 
lution effects in the observations are taken into account, this be¬ 
haviour can be reproduced by simulations with plausible amounts 
of AGN feedback. 

The results of IPlanck Collaborati^ i2013tl and 
lAnderson et al.l ( 1201 Sh are intriguing because they seem to 
indicate the current state of a significant fraction of the baryons 
which had previously been “missing”. It is, however, a matter of 
concern that the mass calibration of their scaling relations rests on 
a specific simulation of the formation and evolution of the galaxy 
population. The goal of this paper is four-fold : (i) to check the 
halo mass calibration provided by the model of IGuo et al.l ( l2013l) 
by comparing its predictions with weak gravitational lensing 
measurements for the LBG sample used in the stacking analyses; 
(ii) to transfer the mass calibration of these scaling relations to one 
coming directly from the lensing measurements, thus substantially 
reducing its model dependence; (iii) to evaluate the residual model 
dependence of this calibration and its origin by comparing results 
for a large number of different galaxy population simulations; 
and (iv) to combine the resulting estimate of residual systematic 
uncertainties with uncertainties coming from noise in the lensing 
measurements to obtain realistic error bars on the parameters of 
the recalibrated scaling relations. 

The plan of our paper is the following. In Sec.[2 we describe 
the observational samples and the methods we use to measure dif¬ 
ferential lensing surface density profiles and two-point autocorre¬ 
lation functions. Our N-body simulations and the galaxy forma¬ 
tion models based on them are introduced in Sec. [3] We present 
results for differential lensing surface density profiles for SDSS 
LBGs in Sec. |4] and compare them with predictions from mod¬ 
els with varying cosmologies, N-body realisations and physical 
models for galaxy formation. Galaxy clustering measurements for 
the SDSS LBGs are compared with various models in Sec. [5] Fi¬ 
nally, we r ecalibrate the scaling re lations of IPlanck Collaborationl 
( l2013h and lAnderson et alj ( 1201 Sh in Sec. providing new rela¬ 
tions which account consistently for the uncertainties both in ef¬ 
fective halo mass and in stacked SZ and X-ray flux for each bin of 
LBG stellar mass. In a companion paper. iMandelbaum et alj ( l2015l) 
compare lensing signals for active and passive (e.g. blue and red) 
subsamples of our LBG sample, finding red LBGs to have sig¬ 
nificantly more massive halos than blue ones of the same stellar 
mass, w hile the same conclusion is reached bv IZu & Mandelbauml 
(l2015bll through halo modelling to both lensing and clustering 
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Figure 1. A scatter plot of our sample of locally brightest galaxies in the 
redshift - stellar mass plane. Red equidensity contours enclose 10%, 30%, 
60% and 90% of the galaxies. To avoid overcrowding only 10% of the sam¬ 
ple galaxies are plotted. 


measurements. When quoting observational results, we adopt as 
our fiducial cosmological mo del the first-year Planck cosmology 
dPlanck Collaborationl 12014a . with present values of the Hubble 
constant Ho = 67.3 km s“^/Mpc, of the matter density Qm ~ 
0.315 and of the cosmological constant Da = 0.685). Our simula¬ 
tions assume a variety of other cosmologies as described in Sec. [3] 
In the following we will define the reduced Hubble parameter h as 
h = Ho/100 km s“^/Mpc. 


2 OBSERVATIONAL SAMPLES AND TECHNIQUES 
2.1 The locally brightest galaxy sample 


In order to define a large sample of galaxies that are almost all cen¬ 
trally located and dominant within their dark matter halos, we se¬ 
lect locally brightest galaxies (hereafter LBGs ) from the New York 
University Value Added Galaxy Catalogue ( Blanton et al.l 1200^ 
VAGC), which is b ased on the seventh da ta release of the Sloan 
Digital Sky Survey dAbazaiian et aD l2009L SDSS/DR7). This par¬ 
ent catalogue contains 602,251 galaxies with high quality spectra 
and is flux limited at r = 17.7 (r-band, extinction-corrected Pet¬ 
rosian magnitude). 

LBGs are selected to be galaxies that are brighter in r than 
all other catalogue members projected within 1.0 Mpc and with 
redshift differing by less than 1000 km/s. The SDSS spectroscopic 
sample is incomplete due to exclusion effects when placing fibres 
in crowded fields. To ensure that no candidates have brighter com¬ 
panions wi thout spectroscopy , we have used a photometric redshift 
catalogue dCunha et alJl200Sl . photoz2) based on SDSS/DR7 pho¬ 
tometry to look for additional companions missed in the spectro¬ 
scopic sample. We eliminate any candidate with a companion in 
this catalogue of equal or brighter r-band magnitude and projected 
within 1.0 Mpc, unless the photometric redshift distribution of the 
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potential companion is inconsistent with the spectroscopic redshift 
of the candidate^ 

After this selection 279,343 LBGs remain. Of these, 161,791 
have stellar mass greater than 10^°'® Mq (see Sec. 12.21 for details 
about how the stellar masses are measured). This is the limit to 
which SZ and X -ray signals have been c learly detec ted in the stack¬ 
ing an alyses of IPlanck CollaborationI i2013h and lAnderson et al.l 
( l2015h . hereafter P13 and A15, respectively. The LBG sample 
used here differs slightly from that used in these earlier stud¬ 
ies, despite being selected according to the same criteria, be- 
cause in this paper we interp ret all data using the cosmology of 
IPlanck Collaborat ionI 1 2014ar) whereas the earlier papers adopted 
that of WMAP7 jKomatsu et al.ll20113) . This slightly changes the 
isolation criterion in projected separation, and also increases the 
stellar mass assigned to each SDSS galaxy by about 10% (because 
of the switch from h = 0.704 to h = 0.673). 

Redshift and colour distributions of the selected LBG sample 
are shown in Fig. 1 of PI3. These are very similar to the corre¬ 
sponding distributions for the parent sample, except that the blue 
fraction of LBGs is significantly enhanced for objects with stel¬ 
lar mass below that of the Milky Way, shifting the population to 
slightly higher redshift. Redshift distributions for the specific stel¬ 
lar mass bins used for SZ and X-ray stacking are shown in Fig. 2 
of A15, and the numbers, mean distances, estimated mean halo 
masses and other properties of these subsamples are given in Ta¬ 
ble 1 of that paper. The mean redshift rises from 2 : 0.07 at 

M* ~ 10^®Mq to 2 : ~ 0.3 at M* ~ lO^^M©. In Fig.[T] we show 
a scatter plot of our sample in the redshift - stellar mass plane. 
Red equidensity contours enclose 10%, 30%, 60% and 90% of the 
points. Clearly more massive LBGs are biased to higher redshifts 
as a result of the SDSS flux limit, and at any given stellar mass 
fewer than half of the objects lie at redshifts where the sample can 
be considered complete. We will discuss in Sec. [231 and Sec. 13.61 
how we account for such incompleteness. 

P13 and A15 use d a galaxy forma tion simulation (denoted 
G13-W7' below) from I Guo et al.l ( 12013h to explore the expected 
relation of the LBGs to the dark halos in which they reside. Fig. 2 
of P13 shows that at every stellar mass, more than 83% of LBGs 
should be the central galaxies of their halos - of the remaining 
“satellites” two thirds are included because they are more luminous 
than the central galaxy of their halo. A scatter plot of main halo 
mass against stellar mass in Fig. 3 of P13 demonstrates that while a 
well-defined relation is predicted, it has substantial scatter, a signif¬ 
icant high-mass tail due to the remaining satellites, and a noticeable 
offset between star-forming and passive systems. Histograms of the 
expected halo mass at each LBG stellar mass are given in Fig. B.l 
of PI3. These are roughly lognormal with a FWHM of about an 
order of magnitude. Histograms of the projected offsets of satellite 
systems from halo centre are shown as a function of stellar mass in 
Fig. C.2 of P13. These are typically a few hundred kpc. 

The distribution of offsets from the true halo centre for satel¬ 
lite LBGs are presented in Fig. C2 of PI3, where it is also shown 


^ Instead of providing a single value of photometric redshift, the photoz2 
catalogue provides the probability distribution over 0 < z < 1.47. We 
consider the photometry of the neighbour to be inconsistent with its being 
a true companion if the probability of being at a redshift less than or equal 
to that of the candidate is less than 0.1, which means the candidate lies in 
the low redshift tail of the photometric redshift probability distribution for 
the companion, and it is thus unlikely to be a true companion. We did not 
consider the case when the candidate lies in the upper tail of the distribution 
because the spectroscopic galaxies mostly have 2 < 0.3. 


that a more strictly selected sample of LBGs that are locally bright¬ 
est within rp =2Mpc (2000km/s along the line-of-sight) gives al¬ 
most the same distribution of offsets. With this more strictly se¬ 
lected sample of LBGs, the total number has decreased by 30% for 
LBGs that are more massive than logjQ M, = 11. However, the 
sample purity of true central galaxies is only slightly increased to 
about 87% at logjQ M* ~ 11. The selection criteria in P13, A15 
and this paper are thus a compromise between sample size and pu¬ 
rity. The robustness of the scaling relations against changes in the 
isolation criteria has been tested in PI3 using this more strictly se¬ 
lected sample, showing that the scaling relation is insensitive to the 
isolation criteria. 


2.2 Stellar mass estimation and a reference galaxy sample 

The stellar masses used in our analysis were all estimated with 
the K-correct software of iBlanton & Roweii l l2007h by fitting stel¬ 
lar population synthesis models to K-corr e cted, five-band pho¬ 
tometry from SDSS assuming a IChabrieil ( l2003l) stellar initial 
mass function. K-correc0 uses more than 400 spectral templates, 
of which most co r respon d to instantaneous burst models from 
iBmzual^Oiariojj 20031) based on the Padova1994 isochrones 
( Fagotto et alj 199^ B B). 

iMoustakas et al.ll2013l) compared these stellar masses to in¬ 
dependently derived values based on 12-band UV to mid-infrared 
photometry. Comparison has also been made with MPA/JHU stel- 
lar masses B based on the HSa and D„(4000) spectral indices 
dKauffmann et al.ll2003l) in addition to SDSS photom etry, as well 
as wit h masses from an extension of this method by I Salim et^ 
( l2007t) . Differences between these methods are almost always less 
than 0.1 dex, with K-correct tending to give systematically slightly 
higher (lower) values than the other schemes at low (high) stel¬ 
lar mass. A fitting formula quantifying the difference between K- 
correc t and MPA/JHU stellar masses is available in iLi & WhTt3 
( I 2 OO 9 I) . Note P13 and A15 used the same K-correct stellar masses 
to derive the effective halo mass for the SZ/X-ray scaling relations. 
Any systematic uncertainties beyond the stellar mass, is supposed 
to be self-consistently removed through the lensing recalibration 
(see Sec. [Hi. 

In Sec. 5 of this paper we analyse the clustering of LBGs by 
computing cross-correlations between stellar mass-limited subsam¬ 
ples and a reference galaxy sample. We choose the latter to be all 
galaxies in the parent spectroscopic sample that are more massive 
than 1 O^°M 0 . 


2.3 The source catalogue for lensing 

The source galaxy catalogue used for our gravitat i onal le nsing 
study has been described in detail in iReves et'H] ( 1201 2h . and 
thus we introduce it only very briefly here. Catalogue constr uc- 
tion uses the re-Gaussianization method dHirata & Seliakll2003l) to 
correct for the point spread function (PSF) which affects the ob- 
served shape of g alaxi es. Thanks to this PSF correctio n scheme, 
iReves et ^ d2012) and iMandelbaum et al.l d2012Ll201^ were able 
to reduce systematic errors significantly compared with previous 


^ http://howdy.physics.nyu.edu/index.php/KcoiTect 

® The data were produced by a collaboration of researchers (currently or 
formerly) from the Max-Planck Institute for Astrophysics (MPA) and the 
John Hopkins University (JHU). 
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SDSS source catalogues. More recently, iHan et af] (1201 .4) devel¬ 
oped an improved maximum-likelihood approach, and with the 
same source catalogue they measured relations between halo mass 
and grou p luminosity with h igh accuracy. The source galaxy cat¬ 
alogue of iReves et alj ( 1201 2h covers 9243 deg^ and fully overlaps 
the SDSS/DR7 footprint within which our LBGs are selected. It 
includes all objects with extinction corrected r-band model mag¬ 
nitudes brighter than 21.8. Photometric redshifts of the sources 
have been estimated using the SDSS hve-band photometry and 
the Zurich Extragala ctic Bayesian Redshift Analyzer (ZEBRA; 
iFeldmann et al.|[2006h . The effect of using photom etric redshifts in 
lensin g analyses like ours has been quantihed by iNakaiima et al.l 

(l2012h . 


2.4 Measuring the differential density profile 

In this study we follow the method desc ribed in iMandelbaum et"^ 
( l2005h and iMandelbaum et all ( l2006h to calculate mean dif¬ 
ferential projected density profiles for the total mass distribu- 
tion stacked around our selected LBGs. We point the reader to 
IMandelbaum et all ( l2005h for full details of this method, and only 
summarise its main features in this section. Some further details are 
given in the companion paper. IMandelbaum et alj ( l2015h . 

The mean differential surface mass density profile, AE(rp), 
is defined as the difference between the mean surface density en¬ 
closed by projected radius rp (denoted E(< r^) and the mean sur¬ 
face density at that radius (denoted E(rp)). The quantity AE(rp) 
can be related to the mean tangential shear, 7 t, and the lensing crit¬ 
ical density, Ec, 


AE = 7tEc, 


( 1 ) 


where Ec is defined as 

Ec = — -(2) 

TttG DisDi 

Di and Da refer to the angular diameter distances of lens and 
source, respectively, and Dis is the angular diameter distance be¬ 
tween lens and source. We use physical separations in our analysis 
rather than comoving separations. 

The mean tangential shear can be related to the directly mea¬ 
surable mean tangential ellipticity, et, of source galaxies, the two 
differing by a factor of twice the shear responsivity, defined as the 
response of the ensembl e averaged ellipticity to a small shear (see 
iBemstein & Jarvisl2002tl . Thus, the mean differential surface mass 
density profile can be estimated from the mean tangential ellipticity 
as follows: 


2i? 


(3) 


where R is the shear responsivity. This is essentially a weighted 
average of etEc over the set of stacking centres. The weight, Wi, 
involves the inverse square of the total noise, composed of mea¬ 
surement noise, Uet.i and shape noise, asN, 


Wi = 




2 

SN 


(4) 


Signals are measured around both real galaxies and random 
points. The signal around random points should be subtracted from 
that around real galaxies to remove contributions from systematic 
shear, though on the scales probed in our analysis such system- 
atics are not significant. Moreover, the signal must be multiplied 
by a factor of B{r) = n(r)/nrand(r), which is the the ratio of 


the number density of sources around real galaxies relative to that 
around random points. This accounts for dilution of the signal by 
physically-associated pairs. In the end, we estimate the errors of 
AE(rp) through jackknife resampling. 

For all the lensing measurements in this paper, we follow the 
convention of this subfield and quote values for rp in physical units; 
this contrasts with the next subsection and with Section 5 where we 
follow the convention of the galaxy clustering subfield and quote 
values for Vp in comoving units. 


2.5 Clustering measurements 


We measure the projected cross-correlation function between 
LBGs a nd our reference galaxy sample using the Landy-Szalay es¬ 
timator dLandv & Szalav fl993h . 




CD{rp,r„) - CRD{rp,r^) - DRc{rp,r^) ^ ^ 

RcRoirp, r-ir) 


(5) 

where CD is the pair count between central LBGs and the refer¬ 
ence sample, CRd is the pair count between LBGs and the random 
sample corresponding to the reference sample, DRc is the pair 
count between the reference sample and the random sample corre¬ 
sponding to LBGs, and RcRd is the pair count between the two 
random samples. Vp and r,r are the projected and the line-of-sight 
components of the (comoving) separations, respectively. 

This estimator thus uses separate random samples correspond¬ 
ing to the LBGs and the reference galaxies. These are constructed 
to have exactly the same sky coverage, number density, redshift and 
magnitude distribution as the corresponding real galaxies, but the 
sky coordinates are randomised within the SDSS footprint. 

The SDSS galaxy samples are flux-limited in the r-band. As a 
result, fainter galaxies are only included at lower redshifts (Fig.ITJ. 
In order to get volume-limited statistics to allow comparison be¬ 
tween samples and with the simulations, we weight individual 
pairs in the above estimator by l/min(Ljnax,i,i, Knax.i.j), where 
Vuia,x,i,i and Vmax,i,j are the total volumes within which galaxies 
i and j could be observed, given the SDSS footprint and appar¬ 
ent magnitude limit. This weighting is adopted for all variations 
of pair counts, i.e. for CD, CRd, DRc and RcRd- It results 
in volume-weighted statistics provided the mean number per unit 
volume of correlated pairs with given physical properties (e.g. in¬ 
dividual galaxy luminosities and stellar masses and pair separation) 
varies at most weakly with redshift over the range 0 < 2 < 0.15 
which contributes most of the weight to the cross-correlation mea¬ 
surements. We have checked that this is indeed the case for the 
models we discuss in the next section. 

After obtaining ^(rp, r^), we integrate it to get the projected 
correlation function. 


Wp{rp) = J ^{rp,r„)dr„, ( 6 ) 

using a depth of ± 40 Mpc/h for the line-of-sight integral. This 
is adequate for projected distances of 10 Mpc or less, as used in 
our analysis (see lPadmanabhan et al.ll200^ . We estimate error bars 
for Wp{rp) through bootstrap resampling 100 subsamples from the 
full LBG sample. Note, for each subsample the contribution to the 
counts from individual LBG-reference galaxy pairs can change due 
to the redistribution of weights assigned to the LBGs. We did not 
bootstrap companions in the reference sample. As noted above, 
when compiling pair counts and quoting results for correlation 
functions we will take rp to refer to a particular comoving projected 
separation. 
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Table 1. Numerical and cosmological parameters for the simulations used 
in our analysis 



MS 

MS-W7 

MS-P' 

box size [Mpc/h] 

500 

500 

480 

softening [kpc/h] 

5 

5 

4.8 

particle mass [Mq/H] 

8.61 X 10® 

9.36 X 10® 

9.62 X 10® 

particle number 

2160® 

2160® 

2160® 


0.25 

0.272 

0.315 

Ha 

0.75 

0.728 

0.685 

fib 

0.045 

0.0455 

0.0487 

Ho [kms“ 1 /Mpc] 

73.0 

70.4 

67.3 

ris 

1 

0.961 

0.96 

0-8 

0.9 

0.807 

0.829 

asnOJ 

0.392 

0.370 

0.415 


3 GALAXY FORMATION SIMULATIONS 

Realistic simulations of cosmic structure formation are very valu¬ 
able when interpreting stacking measurements like those of this pa¬ 
per. To the extent that they represent faithfully both the population 
of objects used as stacking centres and the environments in which 
they form, such simulations can characterise the selection effects 
introduced by observational sample construction and the averaging 
process inherent in the subsequent stacking. P13 and A15 inter¬ 
preted their SZ and X-ray sta cking analyses of LBGs with a sim¬ 
ulation from iGuo et alj ( 1201 3h that used semi-analytic methods to 
follow the formation and evolution of galaxies in a high-resolution 
N-body simulation of a ACDM cosmology. Simulations of this type 
can be carried out in large volumes and tuned to match the proper¬ 
ties of the observed galaxy population quite closely. Here we use 
this and a number of other simulations to assess how the quan¬ 
titative interpretation of our results is affected by uncertainties in 
cosmological parameters and in galaxy modelling assumptions. 


3.1 N-body simulations 

We use two very large dark-matter-only simulations to repre¬ 
sent the dis tribution of dark m atter. These are the Millennium 
Simulation dSoringel et al.l |200^ . MS) which in its original ver¬ 
sion ado pted cosmological p arameters from the first-year data of 
WMAP I Snergel et al. 12003 . WMAPl), and a second simulation 


dGuo et alj2013L MS-W7) which is identical to the MS except that 
the cosmologica l parameters were taken from the seventh year anal¬ 
ysis of WMAP dkomatsu et al.ll20lT5. WMAP7). The fluctuation 
phases for the MS-W7 initial conditions are taken from the pub- 
lic , multi-scale, G aussian, white-noise field (Panphasia) described 
in Ijenkin^ d2013l) . and are different from those used for the MS. 
Table [T] provides basic numerical and cosmological parameters for 
the two simulations. Both follow the same number of particles in a 
cube of the same side length from redshift 127 to the present day. 
The dark matter particle mass is slightly larger in the MS-W7 than 
in the MS. 


3.2 Scaling to other cosmologies 

In addition to the cosmologies adopted for the original MS 
and MS-W7, we would like the first-year Planck cosmol¬ 
ogy dPlanck CollaborationI 120143) to be represented among our 
galaxy formation simulations. The scaling algorithm developed by 


Table 2. Output redshifts used for the MS (columns 1, 3 and 4) and for the 
MS-W7 (column 2) for the vaiious cosmologies adopted in our analysis 


WMAPl 

WMAP7 

WMAP7-scaled 

Planck-scaled 

0 

0 

0 

0 

0.04 

0.04 

0.04 

0.05 

0.09 

0.09 

0.08 

0.08 

0.12 

0.12 

0.12 

0.11 

0.17 

0.17 

0.17 

0.18 

0.24 

0.24 

0.22 

0.22 

0.28 

0.28 

0.27 

0.26 

0.32 

0.32 

0.33 

0.31 


lAngulo & White! ( 1201 Oh and improved bv lAngulo & Hilb^ ( 1201 Sh 
allows the stored output of an N-body simulation to be scaled from 
the cosmology in which it was originally carried out to a differ¬ 
ent cosmology. This procedure involves three steps: rescaling the 
box size, particle mass and velocities, relabelling the output times, 
and rescaling the amplitudes of individual large-scale fluctuation 
modes. For this paper we scale the original MS simulation to both 
the WMAP7 and the Plan ck cosmologies, i n the first case using the 
parameters of Table 1 iniGuo et alj d2013ll. a nd in the second the 
parameters of Fig. 1 of lHenrigues et alf i 2015h . We scale the MS to 
the WMAP7 cosmology to enable direct comparison with MS-W7 
for the properties of interest in this paper. Furt her comparison o f 
the direct and scaled simulations can be found in iGuo et"^ ( l2013h . 

Table lists available simulation redshifts over the range 
needed for our analysis. The outputs we use for the original MS 
and MS-W7 are the same. After scaling the MS to the WMAP7 
and Planck cosmologies, the time coordinate is relabelled and the 
available redshifts are different; for WMAP7, the new 2 = 0 cor¬ 
responds to 2 = 0.28 in the original simulation, while for Planck 
it corresponds to 2 = 0.12. The scaled simulation size and parti¬ 
cle mass, and the new cosmological parameters are given for the 
Planck case in Table[T] For the WMAP7 case, the cosmological pa¬ 
rameters are the same as for MS-W7 but the scaled box size and 
particle mass are 521.55 Mpc/h and 1.0062 x 10® M©/h, respec¬ 
tively. In the next section, we introduce the galaxy formation mod¬ 
els used in our analysis, and all models based on a scaled N-body 
simulation will be indicated by a prime (') after their names. 


3.3 Galaxy formation models 

We have used four recent versions of the Munich semi-analytic 
models to follow the formation and evolution of galaxies in the 
evolving halo/subhalo population of our four N-body simulations 
(the original MS and MS-W7, together with the two scaled ver¬ 
sions of the MS). Altogether, we consider eight different galaxy 
formation simulations, which are summarised in Table [3] Among 
these, Gil - Wl, G13 - W7 and G13 - W 7 ', and H15f - P ' 
are the models pub lished in IGuo et aO (1201 ih . I Guo et aP ( l2013ll 
and lHenrigues et all (l2015h . respectively. For each of these, the un¬ 
certain star formation and feedback efficiencies were tuned to pro¬ 
duce close fits to SDSS data on the stellar mass functions of low- 
redshift galaxies, and the resulting model also reproduces observed 
luminosity and autocorrelation fun ctions quite well. In the case of 
H15f — P'. lHenrioues et all ( 1201 Sl) used the abundances and pas¬ 
sive fractions of galaxies at 0.4 ^ 2 3 as additional constraints 
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Table 3. The galaxy formation simulations used in our analysis, showing their different combinations of astophysical modelling (rows) and background 
cosmology (columns). 



WMAPl 

WMAP7 

Planck 

Guo et al. (2011) 

Gll - Wl 

Gll - W7' 

Gll - P' 

Guo et al. (2013) 


G13 - W7, G13 - W7' 


Henriques et al. (2015)-p 

H15p - Wl 


H15p - P' 

Henriques et al. (2015)-f 



H15f - P' 


on model parameters. Details of how well these models reproduce 
the observed galaxy population can be found in the original papers. 

Our model suite thus includes three different cosmologies, 
WMAPl, WMAP7 and Planck. The Planck cosmology is available 
only by using the scaling algorithm introduced above, whereas for 
the WMAP7 cosmology we have both a direct run and one obtained 
through scaling the MS. These three cosmologies are indicated by 
the columns of Table [3] 

Our model suite includes four different sets of gal a xy for ¬ 
mation paramete rs fo llowing the rec i pes of IGuo et aO ( 1201 ih . 
IGuo et al] ( l2013h and iHenrigues et alj (l2015h . These correspond 
to the rows of Table ID Note that we have two versions of the 
IHenrigues et al.l ( l2015n model. H15p is a preliminary version cor¬ 
responding to that originally submitted and posted on the preprint 
archive, while H15f is the final version revised and retuned in re¬ 
sponse to comments from a referee. Compared with the prelimi¬ 
nary version, the threshold between major and minor mergers was 
shifted to produce a better match to the observed morphology dis¬ 
tribution of galaxies, and a minor problem with the merger trees 
was fixed which had led to unphysically rapid growth of the central 
galaxy of the most massive halo in the simulation. 

Thus, Gil — Wl, Gll — W7' and Gll — P' have exactly 
the same galaxy formation physics and parameters, but different 
cosmologies. G13 — W7 and G13 — W7' are based on the same 
cosmology and galaxy formation physics (including all associated 
parameters) but the former is implemented on MS-W7 and the lat¬ 
ter on a scaled version of the MS. H15p — Wl and H15p — P' 
share the same galaxy formation physics, including all parameters, 
but are based on WMAPl and Planck cosmologies respectively. 
H15f — P' is an updated version of H15p — P', as just noted. 

3.4 Corrections to stellar mass through abundance matching 

Although the published semi-analytic models of the last section 
were all tuned to give reasonable matches to the SDSS stellar mass 
function, a number of other observational constraints such as their 
luminosity functions, their size, morphology and kinematics dis¬ 
tributions, and their clustering were also considered when setting 
up the physical modelling and deter mining paramete r s. In p artic- 
ular, the parameters of the model of IHenrigues et al.1 l l2015h were 
adjusted by MCMC sampling to fit not only the low redshift stel¬ 
lar mass function, but also the abundance and passive fractions of 
galaxies over the redshift range 0 < 2 3. Because of these mul¬ 

tiple constraints, the quality of the fit to the observed low-redshift 
SDSS stellar mass function varies between models. This can signif¬ 
icantly affect the mass of the halos assigned to galaxies of a specific 
stellar mass, particularly at high mass. 

In order to separate differences in model predictions caused 
purely by deviations from the observed stellar mass function, 
from differences caused by cosmology or by other aspects of how 



logM./Mg 


Figure 2. The con'ection to stellar masses, A log M,/Mq, needed to bring 
model stellar mass functions into exact agreement with the observed SDSS 
stellar mass function. The at-axis is the stellar' mass after coiTection, and 
the correction is defined to be negative when it reduces the stellar masses 
assigned to galaxies. A hoi'izontal orange dashed line marks zero coirection 
to guide the eye. 


(sub)halos are populated with galaxies, we use an “abundance 
matching” method to correct the stellar masses in each model so 
that it reproduces the SDSS mass function exactly. We preserve 
the ranking of the model g alaxies in ste l lar ma ss, but shift them to 
match the mass function of iLi & Whitel (l2009l) . Explicitly, we take 
the cumulative number density of galaxies at each value of M* 
from SDSS, and find the stellar mass for which the model predicts 
the same cumulative abundance. The difference between SDSS and 
model stellar masses at this abundance then gives the necessary cor¬ 
rection as a function of stellar mass. Fig. |2] gives the result of this 
exercise. Notice that over most of the mass range shown the cor¬ 
rections are quite small, demonstrating that the galaxy formation 
simulations have indeed reproduced the observations quite accu¬ 
rately. In particular, this is true for G13 — W7^ the model used for 
comparison in the stacking analyses of P13 and A15. 

Halo mass - stellar mass relations before and after these cor¬ 
rections to the stellar mass are shown in the top and bottom panels 
of Fig. respectively. Twelve stellar mass bins have been chosen 
with a width of 0.1 dex in log M*. For each bin, we calculate the 
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Figure 3. Mean halo mass, calculated as (log M 200 ). as a function of stellar 
mass for LBGs in our eight galaxy formation simulations. The upper panel 
shows results from the original simulations, while the lower panel shows 
results after the corrections of Fig. |2] have been used to bring the stellar 
mass function of each model into exact coincidence with that of the SDSS. 

mean log halo mass, (log M 200 ) and plot it against the log mean 
stellar mass, log(M*), for all eight galaxy formation simulation^. 
We choose to use (log M 200 ) instead of log(M 2 oo) throughout the 
paper because, as we will see, halo mass distributions at fixed stellar 
mass are broad and closer to lognormal than to normal. In addition, 
this definition of mean mass is closer to the effective halo mass im¬ 
plied by SZ, X-ray and lensing observations of stacked BCGs than 
is (M 200 ) (see P13, A15 and Sect. below). On the other hand, 
we choose to use log(M*) because for the narrow stellar mass bins 
used in this paper the halo mass distribution changes slowly over 
the bin. We have checked that using (log M*), log(M*) or the mid¬ 
dle bin values gives very similar results. 

For the bottom panel of Fig. [3] the stellar mass of each simu¬ 
lated galaxy has been corrected by an amount taken from Fig. 
The simulations compared here thus have identical stellar mass 
functions. It is encouraging to see that the scatter between simula¬ 
tions, although not entirely eliminated, is significantly less than in 
the top panel. The standard deviation among the eight models in the 
upper panel and for the twelve stellar mass bins are 0.104, 0.090, 

^ For central BCGs, M 200 denotes the mass of a surrounding spherical 
region with mean enclosed density 200 times the critical density of the Uni¬ 
verse. For satellite BCGs, we use the M 200 value of the associated central 
galaxy. 


0.073, 0.060, 0.068, 0.098, 0.120, 0.111, 0.098, 0.109, 0.130 and 
0.131 (from the least to the most massive bin), while for the bottom 
panel the values become 0.065, 0.065, 0.054, 0.039, 0.036, 0.061, 
0.070, 0.056, 0.073, 0.093, 0.136 and 0.198. This demonstrates that 
much of the scatter between the original models can indeed be as¬ 
cribed to imperfect matching of the observed stellar mass function. 
At the most massive end, however, the scatter remains substantial 
and even becomes bigger, while the models split into two groups af¬ 
ter the correction: the Guo et al. group and Henriques et al. group. 
This indicates that while models with similar physical recipes differ 
mainly because of the varying degree to which they fit the observed 
SDSS stellar mass functions, models with differing physical recipes 
can produce different predictions, even when they are applied to the 
same N-body simulation and produce the same stellar mass func¬ 
tion. This is because they are differently sensitive to the detailed 
assembly histories of massive halos. Note Fig. only shows the 
mean halo mass at fixed stellar mass. In Sec. [4] we will show the 
distribution of log M 200 , finding that for models with similar mean 
halo mass at given LBG stellar mass, the shape of the halo mass dis¬ 
tribution is very broad and can differ substantially between models. 

3.5 Selecting locally brightest galaxies in the models 

We select LBGs in our simulations as in P13 and A15 using cri¬ 
teria patterned_closely_on tho se used to select SDSS LBGs (see 
also lWang & Whit3l2012l and IWang et alj[^014h . Specifically, we 
project the simulation box along its z-axis (representing the line- 
of-sight) and assign each model galaxy a redshift based on its z- 
coordinate and velocity. Isolation criteria in projected separation 
and redshift difference can then be applied to the simulation in the 
same way as to the observational data. Such direct projection of the 
simulation box maximises statistical signal, but fails to represent 
effects due to the flux limit of the real survey, to the K-corrections 
needed to obtain rest-frame magnitudes from the observations, to 
the incompleteness for close pairs caused by fibre exclusion ef¬ 
fects, to the complex geometry of the real survey, and to the evo- 
lution of the real pop ulati on across the reds hift range surveyed. 
IWang & Whitel l l2012f) and IWang et al.l (1201 4h have tested that, for 
many of these effects, the results obtained by simply projecting the 
simulation box are unbiased compared to results obtained from a 
full light-cone mock catalogue. As described below, we combine 
data from a range of simulation outputs (see Table |2} to account 
consistently for the effects of evolution and of the flux limit of the 
observational survey. 

In the simulations we know whether each selected LBG is 
truly a central galaxy. As noted in Sect. the fraction of LBGs 
that are central galaxies (i.e. the purity of the sample) is predicted 
by G13 — W7' to exceed 83% at all stellar masses. The maximum 
contamination (-^ 17%) by satellites is predicted to occur at stel¬ 
lar masses slightly above IO^^Mq. P13 checked those LBGs that 
are satellites, finding that for M* ^ lO^^M© about two-thirds are 
brighter than the true central galaxies of their halos. The remainder 
are fainter than their centrals, and are considered locally brightest 
because they are projected more than 1 Mpc from their centrals 
(60%) or have redshifts differing by more than 1,000 km/s. This 
high level of purity allows us to consider our LBG sample as a 
sample of central galaxies with only minor contamination by satel¬ 
lites. 

The projected offsets of satellite LBGs from true halo centre 
are analysed in Appendix C of PI3 and are typically a few hundred 
kpc. It is shown there that a more strictly selected sample of LBGs 
(selected to be locally brightest within a projected separation of 
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Figure 4. Left: Differential surface density profiles centred on LBGs in the Gil — P' model are compared to our observational results for SDSS (filled circles 
with error bars). Solid lines are for the original model, while dashed lines include coiTections to bring the stellar mass function into exact agreement with 
SDSS. Lines and symbols are coloured based on the logarithmic range in stellar mass (in units of Mq) as indicated by the legends. Vp is given in physical (not 
co-moving) units in this and all subsequent plots of this type. The large simulation box size enables us to achieve very small errors for the model (comparable 
to the symbol size). Note we only show the errors for solid lines. Errors for the dashed lines are of comparable size. Right: The host halo mass distributions 
for LBGs at z = 0 within these same stellar mass bins and for the same model, Gil — P'. Vertical lines mark the characteristic halo mass, {log M 200 ) for 
each histogram. Black lines correspond to the original model, and red lines to the model after correction to reproduce the SDSS stellar mass function exactly. 


2 Mpc and a redshift difference of 2000 km/s) has a very similar 
offset distribution and, despite a 30% reduction in sample size, only 
a slightly improved level of purity, about 87% for M* IO^^Mq. 


3.6 Stacked leasing density proiites in the simuiations 

Our SDSS sampie of LBGs is flux-iimited in the r-band and hence 
strongly biased towards higher stellar mass (actually higher lumi¬ 
nosity) at higher redshift (see Fig. [Til. The redshift range spanned is 
approximately 0 < 2 < 0.35, so that high mass LBGs are seen at 
a systematically slightly earlier stage of cosmic evolution, thus on 
a slightly different M* - M 200 relation, than low mass systems. 
Moreover, the lensing signal measured around each galaxy in a 
stack is weighted by the inverse s quare of the (redshif t- depen dent) 
lensing critical density, Ec (see iMandelbaum'et ^ l l200d) and 
Sec. 12.4b . To compare our simulations appropriately to our SDSS 
measurements, we thus need to reproduce the stellar mass - redshift 
distribution of the observations and to assign the same weights to 
objects when stacking. Our procedure to achieve this is as follows. 

(1) Each SDSS LBG is matched to ten simulated LBGs based on 
redshift, stellar mass and luminosity. Explicitly, an SDSS LBG is 
first matched to the simulation output which is closest to its red¬ 
shift. Ten simulated LBGs are then picked from this output, with 


stellar mas and r-band absolute magnitud^fl differing from the 
SDSS values by less than 0.02 dex and 0.05 mag, respectively. If 
fewer than ten simulated galaxies satisfy these requirements, the 
tolerances are increased iteratively by factors of 1.5 until ten sim¬ 
ulated galaxies are matched. In this way, the SDSS and simula¬ 
tion LBG samples end up with identical joint distributions of stellar 
mass, r-band luminosity and redshift. 

(2) Each matched LBG in the simulation is assigned the redshift 
of the corresponding SDSS LBG, and is given a “weight” equal 
to the mean weight (Eqn. II averaged over the background source 
population of SDSS galaxies at that redshift. 

(3) Matched LBGs are grouped by stellar mass in the simulation, 
after scaling by a factor of (hsim.cosm./hpianck)^- This is to treat 
the simulation as though we do not know its true cosmology, and 
so compare it to observation assuming our fiducial Planck cosmol¬ 
ogy. The differential surface density profiles of dark matter parti¬ 
cles around the matched LBGs are calculated as a function of phys¬ 
ical (not comoving) projected radius, Vp, and are averaged using the 
weights assigned in the previous step. 


® In this paper, stellar masses are always quoted in units of the solar mass, 
Mq. The original stellar mass in the VAGC catalogue was given in units 
of M 0 h“^. To eliminate the Hubble constant, we adopt the Planck cos¬ 
mology. The cosmologies of some of our simulations differ, however, so 
we scale the stellar masses (in units of Mq) of simulated galaxies by 
(hsim.cosm./hplanck)^ when matching to the SDSS galaxies. 

^ The term —5 log h in the SDSS magnitudes is eliminated assuming the 
Planck cosmology. For stellar-mass-corrected models, the r-band luminos¬ 
ity is corrected by the same factor as the stellar mass in order to keep the 
stellar mass-to-light ratio unchanged. 
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4 STACKED LENSING SUREACE DENSITY PROFILES 
4.1 Fit quality for a “good” model 

We start by presenting our measurements of differential surface 
density profiles for stacks of SDSS locally brightest galaxies. We 
compare these with predictions from Gil — P', which turns out to 
be the simulation that gives the best overall fit to the observational 
data. Gil — P' is also special in that it requires the smallest cor¬ 
rection of all our models (see Fig. IS to bring its high-mass stellar 
mass function (logM*/M 0 > 10.7) into exac t agreement with 
that derived for the SDSS hv iLi & Whitel l l2009ll . Its behaviour at 
high mass in the lower panel of Fig.[3]is notable for dropping below 
that of the other Guo et al. models. 

Points with error bars in the left-hand plot of Fig. |4] show ob¬ 
servational results from SDSS for LBGs stacked into seven bins of 
log M,/M 0 with boundaries as indicated by the legend. Here and 
in all later similar plots, AE is given as a function of projected ra¬ 
dius, Tp, in physical (rather than comoving) units. Solid lines are the 
predictions for AE(rp) for these same LBG stellar mass bins from 
Gil — P'. Dashed lines are based on the same model, but after cor¬ 
recting to bring its stellar mass function into exact agreement with 
SDSS (see Sec. 13.4b . Because of the large number of dark matter 
particles in the simulation and the large number of LBGs stacked, 
the statistical uncertainties on the model predictions are very small, 
and are indicated only for the solid lines. For SDSS also the un¬ 
certainties are small at most radii in all but the lowest stellar mass 
bins. 

The agreement of Gil — P' with the SDSS measurements is 
very good, both in amplitude and in shape over almost three orders 
of magnitude in projected radius and 1.5 orders of magnitude in 
stellar mass. Both theory and observation show a clear transition 
between a small-scale “one-halo” term, where the signal is domi¬ 
nated by the LBGs’ own halos, and a large-scale “two-halo” term, 
where the signal is dominated by the correlated environment sur¬ 
rounding the LBG halos. The transition region is distorted by a 
feature at ~ 1 Mpc induced by our requirement that LBGs should 
have no brighter companion projected within this radius. The stellar 
mass corrected and uncorrected cases are very similar at high mass, 
where the corrections are small; differences are more noticeable in 
the lowest mass bins. For the 18 data points shown, the differ¬ 
ences between the observed profiles and the uncorrected model are 
20.45, 52.40, 14.95, 36.48, 44.79, 48.93 and 54.79 for the most 
massive to the least massive bin. For the stellar mass corrected 
model, the corresponding values are 32.74, 44.14, 12.29, 36.10, 
54.03, 67.35 and 66.12. While many of these values are formally 
sufficient to exclude the model, they actually represent remarkably 
good fits given the small size of many of the formal (bootstrap) er¬ 
ror bars and the fact that no model parameters were adjusted to fit 
these data. 

As we will see later through comparison with other simula¬ 
tions, the excellent agreement of Gil — P' with SDSS is especially 
notable for the three most massive bins. The model also agrees rea¬ 
sonably well with the observations for the four lower mass bins, 
but it seems to have some difficulty reproducing the one-halo/two- 
halo transition and the 1.0 Mpc/h feature induced by our isolation 
criterion. The one-halo/two halo transition has often proved diffi¬ 
cult to fit accurately with physical or phenomenological models of 
clustering (see, e.g .jZehav i et al.|[20ll] ; Ivan den Bosch et ^l2013l ; 

IZu & Mandelbaumll2015aL for Halo Occupation Distribution fits to 

galaxy autocorrelations in SDSS). The 1.0 Mpc/h feature is due to 
the requirement of no brighter companions within 1 Mpc, which 
translates to a plummeting of the differential surface density pro¬ 


files at this scale because galaxies trace the underlying mass distri¬ 
bution. This is especially true for low mass LBGs, as our method 
makes us preferentially select rather isolated systems. 

The right-hand plot in Fig. |4] shows histograms of host halo 
mass for LBGs in our seven stellar mass bins. For simplicity, here 
and throughout the paper we use LBGs selected in the z = 0 snap¬ 
shot to make these plots. For non-central LBGs, the M 200 val¬ 
ues adopted are those of the central galaxies of which they are 
satellites. Black and red histograms in Fig. |4] refer to the stellar 
mass uncorrected and stellar mass corrected cases, respectively, as 
noted in the legend. Vertical lines mark the characteristic halo mass, 
(log M 200 ), for each bin. We adopt (log M 200 ) (instead of, for ex¬ 
ample, log(M 2 oo)) because the halo mass distributions are broad 
and roughly lognormal so that (logM 2 oo) values are in all cases 
close to the central values of the distributions. They also turn out to 
be close to the “effective” halo masses defined by the SZ and X-ray 
stacking analyses of PI3 and A15. 

The scatter in halo mass at fixed stellar mass is large and de¬ 
pends both in width and in shape on stellar mass (and, as we will see 
below, on galaxy formation model). For the four most massive bins, 
the corrected and uncorrected models give almost the same distri¬ 
bution, explaining why the solid and dashed lines are almost in¬ 
distinguishable in the left plot. For the three least massive bins, the 
distributions shift slightly to the right for the mass corrected model, 
giving higher characteristic halo masses and higher predictions in 
the left-hand plot. The amplitude differences in the predicted tens¬ 
ing signal agree qualitatively with the shifts in characteristic mass 
of the halo mass distributions. Despite these small amplitude shifts, 
both solid and dashed curves agree quite well with SDSS. 

4.2 Fit quality for two published Munich models 

The model analysed in the last section, Gil — P', impl e ments the 
galaxy formation physics and parameters of IGuo et alj ( 1201 ih on 
a version of the MS sc aled to the Planck cosmology first used in 
iHenrigues et al.l ( l2015h . As a result, it does not correspond to any 
previously published galaxy formation simulation. In order to start 
our exploration of the model-dependence of predicted lensing sur¬ 
face density profiles, we now present results for two published and 
publicly available simulations, Gll — Wl, which is based o n the 
WMAPl cosmology and was published bv iGuo et alj 1 I2OI ih . and 
H15f — P', which is based on th e Planck cosmology and was pub¬ 
lished by IHenrigues et all i2015tl . Fig.|5]shows differential lensing 
surface density profiles and halo mass distributions for these two 
models in exactly the same format as Fig.|4] Both stellar mass cor¬ 
rected and stellar mass uncorrected results are shown for both mod¬ 
els. Note round dots with errors are lensing profiles based on SDSS, 
which are exactly the same in all relevant plots from Fig.|4]to Fig.[ 8 ] 
and also for Fig. [TO] For these figures we always only show errors 
for solid curves, while errors associated to dashed curves have sim¬ 
ilar size. In most cases the model predictions have extremely small 
errors comparable to the symbol size. 

The left side of Fig.[5]shows results for Gll — Wl. As can be 
seen from Fig .|2 this model requires significant and negative cor¬ 
rections at high mass to bring its stellar mass function into agree¬ 
ment with SDSS, i.e. it overpredicts th e abundance of high mass 
galaxies (see Fig. 7 of IG uo et alj|201 ih . After correction, the pre¬ 
dicted surface density profile for the two highest mass bins lies sig¬ 
nificantly above the observations, showing that the excellent agree¬ 
ment found for the original model is, in fact, a coincidence. For the 
lower stellar mass bins, the solid and dashed curves are very simi¬ 
lar since the required mass corrections are very small in this stellar 
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logM2(|„[Mg/h] logMjoolMg/h] 

Figure 5. Top Left and Bottom Left: As Fig.|4]but for Gil — Wl. Top Right and Bottom Right: As Fig.|4]but for H15f — P'. 


mass range. For all seven bins, the differences in characteristic halo 
mass (the shifts between the red and black vertical lines in the lower 
plot) agree quite well with the shifts in amplitude between solid and 
dashed lines in the upper plot. This shows that correcting the model 
stellar masses induces a (logarithmic) shift in the “effective halo 
mass” (defined as the halo mass which would produce the same 
signal when all LBGs are assumed to have identical halos) which is 
very similar to that in the characteristic halo mass (log M 200 ). This 
implies that the shape of the halo mass distributions at fixed stellar 
mass is not much affected by the correction. 

There is quite good agreement with SDSS for the five lower 
mass bins, except that magenta curves underpredict the SDSS sig¬ 
nal at projected radii of a few hundred kpc. values for the fit 


of the uncorrected Gil — Wl model to the SDSS data are 18.05, 
17.46, 22.85, 49.02, 24.86, 32.48 and 44.18 from most to least mas¬ 
sive bin. For the mass-corrected model, the corresponding values 
are 157.40, 98.22, 14.39, 46.55, 23.26, 35.56 and 47.75. It is clear 
that, after stellar mass correction, the Gll — Wl model fits obser¬ 
vation much less well than Gll — P', at least for the two highest 
stellar mass bins. 

The right side of Fig.^shows similar results for H15f — P'. 
The SDSS signal is clearly underestimated between 200 kpc and 
2 Mpc in the most massive bin, and this problem is not affected 
by the mass correction because this particular model fits the high- 
mass tail of the SDSS stellar mass function very well (see Fig. 2 of 
iHenrigues et al.l2015L note also that the abundance at lower masses 
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Figure 6. Top Left and Bottom Left: As Fig.|4]but now compaiing Gil — W1 and H15p — W1 after both have been corrected to reproduce the SDSS stellar' 
mass function. These two models thus have exactly the same stellar mass function and cosmology (WMAPl) and ai'e implemented on the same simulation 
(the original MS), yet they make substantially different predictions for the lensing profiles and have different halo mass distributions at given stellar mass. Top 
Right and Bottom Right: As the left-hand plots, but now for Gll — P' and H15f — P', two simulations based on the Planck cosmology. 


is underpredicted, giving rise to the positive corrections seen for 
this model in Fig. [2]l. For the five lowest stellar mass bins, sig¬ 
nificant stellar mass corrections cause both the characteristic halo 
masses and the shape of the halo mass distributions to shift. As 
a result, the lensing profile predictions in the upper plot are low¬ 
ered and come into better agreement with SDSS. The values are 
25.34, 25.64, 69.64, 94.26, 165.07, 96.06 and 78.50 for the mass- 
uncorrected model and 25.49, 28.66, 18.17, 55.73, 78.25, 57.83 and 
58.25 after mass correction (again from most to least massive bin). 
Once more it is true that the changes in amplitude of the surface 


density profiles (i.e. the shifts in the effective halo mass for lens¬ 
ing) agree well with the shifts in characteristic halo mass seen in 
the lower plot. 

It is clear that overall Gll — W1 and H15f — P' agree less 
well with SDSS than Gll — P', especially for the few most 
massive bins and after correction to bring agreement with the 
SDSS stellar mass function. For lower stellar mass, the abundance¬ 
matching corrections improve the agreement with observation, at 
least for these models. 
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Figure 7. As Fig.|4] but now comparing G13 — W7 and G13 — W7', two models with the same galaxy formation physics and cosmology {WMAP7), one 
implemented on a simulation executed with that cosmology (MS-W7) and the other implemented on a scaled version of the MS. CoiTections to bring the 
stellar mass function into exact agreement with SDSS have been included in both cases. Differences here are due to cosmic variance and to inadequacies of 
the scaling procedure. 
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Figure 8. As Fig.|4] but now comparing three simulations with the same galaxy formation physics and parameters, and implemented on the same N-body 
simulation (the MS) but with different cosmologies (Gil — P', Gil — W1 and Gll — W7')- Corrections to bring the stellar mass function into exact 
agreement with SDSS have been included for all three models. Differences here are thus due almost entirely to the differing cosmologies. 


4.3 Models differing only in galaxy formation physics 

Fig. [^compares differential surface density profiles and halo mass 
distributions for two galaxy formation models implemented on the 
original MS (Gll — W1 and H15p — W1 on the left side) and for 
two implemented on the MS after scaling to the Planck cosmol¬ 


ogy (Gll — P' and H15f — P' on the right side). Thus the models 
in each pair have the same evolving dark matter distribution but 
assume different galaxy formation physics. All four models have 
been corrected to bring their stellar mass distributions into exact 
agreement with SDSS and so have identical stellar mass functions. 
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Differences within each pair are thus due solely to differences in 
how galaxies of a given stellar mass are assigned to (sub)halos. 

A comparison of Gil — W1 and H15p — W1 shows that de¬ 
spite their identical mass distributions and (corrected) stellar mass 
functions, their surface density profiles differ substantially in the 
most massive bins, Gll — W1 overpredicting the observed signal 
and H15p — W1 underpredicting it. Corresponding values are 
73.34 and 43.42 for the two most massive bins in H15p — W1 and 
157.40 and 98.22 for these same bins in Gll — Wl. The val¬ 
ues for the five less massive bins of H15p — Wl are 35.44, 35.13, 
29.02, 29.64 and 45.66, which are closer to the corresponding 
values for Gll — Wl (see Sec. l4.2t . 

These substantial differences highlight the importance of scat¬ 
ter in the stellar mass - halo mass relation. This scatter is not, of 
course, purely stochastic, but rather a consequence of the scatter 
in assembly history among halos of given mass, which impacts the 
properties of central galaxies differently in different galaxy forma¬ 
tion models. Scatter effects are particularly important at high mass 
where the dependence of cen t ral galaxy prope r ties on halo mass is 
weak (e.g. iReves et al.ll200^ : I Guo et allboid : iMoster et al.ll2010l : 
IWang& Jinall2010h . As a result, central galaxy stellar mass is a 
relatively poor and model-dependent predictor of halo mass on the 
scale of rich clusters. More precise mass proxies, perh aps, cluster 
richness, total luminosity. X-ray or SZ properties (e.g. IRozo et ^ 
l201(]l : |pranck Collaborationll2014bl : [oid et aLlJl2015h will be needed 
to derive robust and competitive cosmological constraints from the 
cluster population. 

The differences between models are very clear in the halo 
mass distributions in the lower panels of Fig. In a given stellar 
mass bin, the distributions can differ substantially in shape, width 
and characteristic mass, with effects being particularly large in the 
higher stellar mass bins. Gll — Wlpredicts less scatter at the mas¬ 
sive end of the halo mass - stellar mass relation than H15p — Wl. 
In the lower stellar mass bins, the characteristic halo masses agree 
better, but the shape and amplitude of the halo mass distributions 
still differ noticeably. 

Similar but less dramatic differences are seen for the two mod¬ 
els based on the MS scaled to the Planck cosmology. In the two 
most massive bins, Gll — P' predicts a significantly stronger tens¬ 
ing signal than H15f — P', although the difference is not as big as 
between Gll — Wl and H15p — Wl. The values for the two 
most massive bins are 25.49 and 28.66 for H15f — P' and 32.74 
and 44.14 for Gll — P'. The behaviour for the five less massive 
bins parallels that for the WMAPl models, with x^ values of 18.17, 
55.73, 78.25, 57.82 and 58.25 for H15f - P'. Although the charac¬ 
teristic halo masses and differential surface density profiles differ 
very little between the two models at these stellar masses, the shape 
and amplitude of the halo mass distributions again differ quite no¬ 
ticeably. 

Despite the substantial shape and and scatter variations be¬ 
tween these models, the shifts in normalisation of the differential 
surface density profiles correspond moderately well to the shifts in 
characteristic halo mass, (logM 2 oo} for all stellar mass bins and 
in both cosmologies. The relative amplitude of the lensing profiles 
provides an estimate of the change in effective halo mass between 
two models (see Sec. for details of how we estimate this). The 
logarithmic shifts from H15p — Wl to Gll — Wl are estimated 
to be 0.469, 0.313, 0.074, -0.076, -0.058, -0.010 and 0.019. from 
the most to the least massive stellar mass bin. For comparison, the 
logarithmic differences between the characteristic halo masses in¬ 
dicated by the vertical lines in the bottom left plot are 0.256, 0.216, 
0.140, -0.023, -0.052, 0.024 and 0.065. Similarly, the shifts ineffec¬ 


tive halo mass from H15f — P' to Gll — P' estimated from their 
differential lensing surface density profiles are 0.238, 0.111, 0.027, 
-0.014, -0.025, 0.001 and 0.028, while the corresponding logarith¬ 
mic shifts in characteristic mass between their halo mass distribu¬ 
tions are 0.150, 0.104, 0.065, -0.012, -0.042, 0.012 and 0.028. The 
shifts in the effective halo mass for lensing and in (log M 200 ) are 
similar but certainly not identical. The detailed shape of the halo 
mass distribution at fixed stellar mass thus affects the ratio of these 
two masses, rendering it model-dependent. 

4.4 True versus scaled N-body simulations 

In this subsection we compare G13 — W7 and G13 — W7'. These 
two simulations have identical galaxy formation physics, includ¬ 
ing all parameters, and assume the same cosmology. However, 
G13 — W7 is based on the MS-W7, which was set up and carried 
out assuming a WMAP7 cosmology, whereas G13 — W7' is based 
on the MS (carried out assuming WMAPl parameters) after rescal¬ 
ing to a WMAP7 cosmology. We can already see from Fig.|2]that 
the two models do not produce exactly the same stellar mass func¬ 
tion; their mass correcti on curves differ b y about 0.05 dex on all 
scales (see also Fig. 5 of iGuo et alj|2013h . Such differences could 
arise from cosmic variance (since the two simulations have inde¬ 
pendent initial fluctuation fields), from the (small) differences in 
time and mass resolution between MS-W7 and the scaled MS, and 
from inadequate accuracy in the scaling algorithm. 

Fig. |7] presents a direct comparison between G13 — W7 and 
G13 — W7' in the same format as Fig. |4] For both simulations 
we have included the corrections needed to bring their stellar mass 
functions into exact agreement with SDSS. With this source of dif¬ 
ference removed, they predict very similar differential surface den¬ 
sity profiles. For the few most massive bins, G13 — W7 (dashed) is 
slightly higher than G13 — W7' (solid) on small scales. The differ¬ 
ence shows up in the x^ values for the two models: 77.42, 38.36, 
13.64, 35.69, 53.71, 47.40 and 50.98 for G13 - W7 and 47.02, 
21.92, 23.69, 38.77, 51.74, 51.59 and 54.48 for G13 - W7'. Thus 
G13 — W7' agrees better with SDSS for the two most massive 
bins, while for the remaining bins the two models have similar 
values. We note that G13 — W7' is the simulation that was used by 
P13 and A15 to interpret their SZ and X-ray stacking analyses, and 
that its agreement with the SDSS lensing data is very good, only 
slightly worse than that of Gll — P', the “best” model of Sect. 14. II 

The ha lo mass functions of the two simulations are compared 
in Fig. 1 of IGuo et al.l ( l2013h . where MS-W7 can be seen to have 
slightly higher halo abundances at z = 0 than the scaled MS. On 
the other hand, Fig.j^shows that at a given stellar mass, G13 — W7 
predicts slightly lower galaxy abundances than G13 — W7'. Thus 
the halo mass - stellar mass relation shifts between the two mod¬ 
els, with halos of a given mass forming central galaxies of slightly 
lower stellar mass in the scaled simulation, perhaps as a conse¬ 
quence of its slightly more widely spaced sequence of output times. 
Correcting both simulations so that they reproduce the SDSS stellar 
mass function removes much of this difference. 

Halo mass distributions are shown on the right-hand side of 
Fig.lH For every stellar mass bin, the two models have similar char¬ 
acteristic halo masses, even though in some cases the shapes of the 
distributions differ significantly. Here again the logarithmic shifts 
in characteristic halo mass are similar in size to those in the effec¬ 
tive mass corresponding to the lensing profiles. From G13 — W7' 
to G13 - W7 the latter shift by -0.063, -0.080, -0.072, -0.053, 
0.003, -0.028 and -0.036, while the former shift by -0.040, 0.055, 
0.039, 0.012, -0.008, -0.020 and -0.032. The two sets of numbers 
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Table 4. Leasing calibrated characteristic halo masses and their systematic uncertainties as inferred using our full set of galaxy formation simulations. We 
use each simulation both in its original form and after correction to bring its stellar mass function into exact agreement with the SDSS, resulting in a total of 
sixteen simulations. For each simulation and each stellar' mass bin, we calculate the mean and median values of log M2oo- These are then calibrated by shifting 
logarithmically by the amount needed to bring the simulated differential surface density profile into agreement with the SDSS data (see Sec.[ 6 ). Column 1 lists 
the boundaries of the stellar mass bins; columns 2 and 4 provide lensing calibrated estimates of (log M 200 ) and median(log M 200 ) for our SDSS sample, 
averaged over all 16 models, while columns 3 and 5 give the rms scatter among estimates for these quantities based on the individual models. M 200 is in units 
ofM 0 h“^ throughout. 


logM*/M 0 

((logM 20 o}) 

o-((logM 2 oo}) 

(median(log M 200 )} 

cr(median(log M 200 )) 

10.-10.4 

11.837 

0.027 

11.778 

0.026 

10.4-10.7 

11.836 

0.022 

11.743 

0.017 

10.7-11. 

12.234 

0.054 

12.131 

0.051 

11.-11.2 

12.709 

0.043 

12.620 

0.056 

11.2-11.4 

13.015 

0.047 

12.975 

0.060 

11.4-11.6 

13.475 

0.064 

13.471 

0.060 

11.6-15. 

13.996 

0.083 

14.000 

0.099 


do not correspond exactly, showing once more a dependence of the 
ratio of effective to characteristic halo mass on the detailed shape 
of the halo mass distributions. 


4.5 Variations with cosmology 

Fig-Elshows differential surface density profiles and halo mass dis¬ 
tributions for three stellar-mass-corrected models with the same 
galaxy formation physics (including all parameters) but different 
cosmologies, Gil — P', Gll — W1 and Gll — W7'. The differ¬ 
ential density profiles for these models are inconsistent with each 
other, showing that cosmology also has a significant effect on the 
lensing predictions. However, understanding the source of these 
differences is not straightforward. Given the extremely small er¬ 
rors of the model predictions (comparable to the symbol size), 
Gll — W1 predicts significantly stronger lensing signals than the 
other two models for the two highest stellar mass bins. For most of 
the points over the whole radius range, the significance is a few tens 
of a with respect to the model errors. Although it is tempting to as¬ 
cribe this to its higher value of ag, Gll — P' has a higher value 
of Got. and simple considerations predict the abundance of mas¬ 
sive halos at « = 0 to depend on the combination which 

is largest for Gll — P' and smallest for Gll — W7' (see Table|T]l. 
At lower stellar masses, Gll — W1 and Gll — W7' both predict 
weaker lensing signals than Gll — P'. These mass-dependent dif¬ 
ferences hinder any simple interpretation of the cosmology depen¬ 
dence of the lensing signals. 

Despite these complications, we again see a fair correspon¬ 
dence between the differences in effective halo mass (inferred from 
the relative amplitude of differential surface density profiles) and 
the differences in characteristic halo mass (found directly from the 
halo mass distributions). The logarithmic shifts in effective halo 
mass from Gll — W1 to Gll — P' are found from their lensing 
profiles to be -0.156, -0.116, 0.069, 0.135, 0.122, 0.101 and 0.094. 
The corresponding shifts in characteristic halo mass are -0.081, - 
0.056, -0.014, 0.074, 0.099, 0.104 and 0.100. From Gll - W7' to 
Gll — P', the shifts in effective halo mass are-0.038,0.029, 0.123, 
0.160, 0.165, 0.151 and 0.164, while the shifts in characteristic halo 
mass are 0.013, 0.045, 0.070, 0.110, 0.106, 0.099 and 0.102. Again, 
the failure of these sets of numbers to correspond perfectly indi¬ 
cates a model-dependence of the ratio between the effective halo 
mass for lensing and our characteristic halo mass, (log M 2 oo)- 


4.6 The model dependence of effective halo masses 

The differential surface density profiles that we have estimated for 
the SDSS through gravitational lensing and calculated directly from 
our galaxy formation simulations are a weighted average of those 
of the individual galaxies in each stellar mass bin. Similarly, the 
stacked SZ and X-ray signals estimated from Planck and ROSAT 
data by P13 and A15 are weighted averages of the signals from in¬ 
dividual galaxies, which the modelling relates to halo mass through 
power-law scaling relations. In all three cases, the stacked signal is 
very similar to that which would be found if all galaxies had halos 
of one particular mass, the effective halo mass for that stellar mass 
bin. The effective halo masses will not be the same in the three 
cases, however, because of differing weights and differing sensitiv¬ 
ity of the observed signal to redshift and halo mass. 

If the distribution of halo masses at given stellar mass were 
always the same shape, for example, lognormal with known width, 
then the ratios of the effective halo masses to each other and to other 
characteristic masses such as {M 200 ) or median(M 2 oo) would al¬ 
ways be the same. The direct measurement of any one of these 
quantities could then be used to calibrate the others. We have seen 
in previous sections, however, that this is not the case; the shape of 
the halo mass distribution varies between stellar mass bins and from 
model to model. This introduces model-dependence into any grav¬ 
itational lensing based recalibration of the scaling relations linking 
the SZ and X-ray properties of halos to their mass, a primary goal 
of this paper. In this subsection we use our suite of galaxy forma¬ 
tion simulations to estimate the level of uncertainty resulting from 
this model-dependence. 

In their Table B.l, P13 give SZ effective halo masses together 
with mean and median values of M 200 for each stellar mass bin 
used in their SZ stacking analysis. A15 give X-ray effective halo 
masses for these same bins in their Table A. 1. In both cases, the ef¬ 
fective halo masses are estimated using G13 — W7' and are much 
closer to median(M 2 oo) than to {M 200 ) in each bin. Since we do 
not wish to repeat the detailed modelling of the SZ or X-ray anal¬ 
yses in this paper, we will assume that the model dependence of 
the ratio of SZ (or X-ray) effective mass to gravitational lensing 
effective mass is similar to that of the ratio of median (M 200 ) to 
gravitational lensing effective mass. The rms scatter in (the loga¬ 
rithm of) this ratio over our set of galaxy formation simulations 
then provides an estimate of the systematic uncertainty induced by 
model-dependence when recalibrating the scaling relations of PI3 
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and A15 based on a single simulation, G13 — W7^ Note that the 
accuracy of this estimate is not explicitly tested in our analysis. It 
can only be tested in detail by running the full observational mod¬ 
elling analysis of P13 and A15 on all the other models, which is not 
possible within the author collaboration of this paper. 

In this paper, we consider eight different galaxy formation 
simulations (see Table[3j, for each of which we create mock SDSS 
samples using both the stellar masses originally produced by the 
simulations and stellar masses corrected so that the simulated stel¬ 
lar mass function matches the SDSS function exactly. Thus, for 
each of our seven stellar mass bins we end up with 16 differ¬ 
ent halo mass distributions with corresponding lensing signal pre¬ 
dictions. For each case and each bin, we find the shift in halo 
mass, A log M 200 , needed to bring the predicted lensing signal into 
agreement with our SDSS measurement (see Section 6 for details) 
and we apply this shift to median (log M 200 ) and (logM 2 oo)- If 
the halo mass distributions predicted for the bin by the 16 different 
models were all the same shape, the 16 values obtained for each 
characteristic mass would coincide. This is not the case, showing 
that the ratio of lensing effective mass to characteristic mass is vary¬ 
ing from model to model. 

Table|4]gives the results of this exercise. For each stellar mass 
bin, column 2 gives the mean over our 16 simulations of the recal¬ 
ibrated value of (logM 2 oo), while column 3 shows the rms scat¬ 
ter among these values. Columns 4 and 5 are similar but for the 
recalibrated values of median (log M 200 ). These means are thus 
our best lensing-based estimates of these particular characteristic 
halo masses for real SDSS galaxies. The scatter values give es¬ 
timates of the systematic uncertainty in these means due to their 
model-dependence. For both quantities the model-dependence is 
quite small (less than about 0.1 dex) but is an increasing function 
of stellar mass, reflecting the fact that central galaxy stellar mass 
becomes a poorer proxy for halo mass in more massive systems. As 
remarked above, we will later adopt the scatter values given here as 
estimates of the systematic uncertainty due to model-dependence 
when using our gravitational lensing observations to recalibrate SZ 
and X-ray scaling relations for our LBGs. 


5 GALAXY CLUSTERING 

In previous subsections we have compared our SDSS measure¬ 
ments of differential surface density profiles around LBGs to pre¬ 
dictions from a range of galaxy formation simulations, finding the 
latter to be significantly affected by the assumed galaxy formation 
physics and by the underlying cosmology, even after the model stel¬ 
lar mass functions have been corrected to bring them into exact 
agreement with SDSS. Such effects are particularly important for 
massive objects and reflect the considerable (and model-dependent) 
scatter in the stellar mass - halo mass relation. In this subsection we 
consider clustering of other galaxies around our LBGs and assess 
whether this can provide additional information to help choose be¬ 
tween models. 

We measure the projected cross-correlation, Wp{rp), between 
LBGs and a reference sample, chosen to be all galaxies more mas¬ 
sive than IO^'^Mq in the SDSS, and also in the simulations. Cross¬ 
correlation with a dense reference sample enables us to get better 
signal than for an autocorrelation, and the low mass threshold re¬ 
duces the strength of the feature introduced by our isolation crite¬ 
rion. For SDSS, the correlation is measured as outlined in Sec. [23] 
The weighting results in a relatively low effective redshift for the 
measured signal, with 50% coming from 2 : < 0.08 and 90% from 



Figure 9. Projected cross correlation, Wp(rp), between LBGs and a refer¬ 
ence galaxy sample with logM*/MQ > 10. Solid and dashed lines are 
for Gll — P' and H15f — P', respectively. We only show errors for solid 
lines. These errors are comparable to the symbol size. In both cases stellar 
masses have been corrected to bring the model stellar mass functions into 
exact agreement with SDSS. Solid dots with eiTor bai's are measurements 
of Wp{rp) for real galaxies based on SDSS (see Sect. [23] . Lensing profiles 
and halo mass distributions for these two models can be found on the right- 
hand side of Fig. |6] Note that in this plot Vp denotes comoving (rather than 
physical) separation. 

2 < 0.15 even for our highest LBG stellar mass bin. This reduces 
evolutionary effects on the signal at the price of increasing the cos¬ 
mic variance in the measurement. Corresponding cross-correlation 
functions for the simulations can be measured directly by count¬ 
ing reference galaxies around LBGs in a projection of the z = 0 
snapshot. Note that, in order to be consistent with convention in 
this subfield, we use rp to denote comoving (rather than physical) 
separation in this subsection. 

In Fig. 13 we compare our SDSS clustering measurement with 
predictions from two representative models, Gll — P' (solid) and 
H15f — P' (dashed), both with stellar masses corrected to bring 
their stellar mass function into agreement with that in the SDSS. 
The lensing signal predicted by these same models is compared 
with our SDSS results in the upper right panel of Fig.|3 The rela¬ 
tive behaviour of the two models is very similar for the two mea¬ 
surements. Gll — P' predicts significantly stronger signals than 
H15f — P' on all scales for the two highest stellar mass bins, given 
the extremely small errors of the model predictions. Note we only 
show errors for the solid lines and these errors are mostly compara¬ 
ble to the symbol size. At lower stellar mass their signals are very 
similar, both for clustering and for lensing. The large-scale clus¬ 
tering amplitudes (rp ^ 2 Mpc/h) measure the biasing of LBGs 
relative to reference galaxies and hence are primarily sensitive to 
their halo mass, which is measured directly on smaller scales by 
the one-halo lensing signal. Fig. |3 shows that Gll — P' predicts 
effective halo masses for the two highest stellar mass bins which 
are slightly above those observed, while those for H15f — P' are 
clearly too low. The large-scale clustering amplitudes of the two 
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models in Fig. [^mirror this behaviour relative to the observations 
as expected. For lower LEG stellar masses, both models fit the 
observed one-halo lensing and large-scale clustering signals quite 
well, with the remaining differences probably within the cosmic 
variance uncertainties of the clustering measurements. 

On smaller scales and for intermediate LEG stellar masses 
(10.7 < logM*/M 0 < 11.4) the observed clustering signal in 
SDSS is significantly larger than in either model, indicating that 
the number of satellites with logM*/MQ > 10.0 is significantly 
underestimated (although, curiously, not around more or less mas¬ 
sive LEGs). We have explicitly checked all the other models, and 
found such small scale excess persists across cosmologies. This is 
unexpected, both because Fig. indicates that the models repro¬ 
duce the mass distribution arou nd these LEGs quite ac curately on 
these same scales, and because IWang & Whitel l l2012h found that 
the count of satellites more massive than log M./M© = 10.0 pro¬ 
jected within 300 kpc of a large sample of bright isolated galaxies 
with 10.8 < logM*/MQ < 11.4 is reproduced quite precisely 
by Gil — W1 (see their Fig. 6 ). We have not isolated the origin 
of this discrepancy, which may in part reflect the underestimation 
of cosmic variance uncertainties by our bootstrap error estimation 
procedure, but it appears unlikely to affect the main topics of this 
paper. Note the apparent turnover on very small scales and espe¬ 
cially for the few most massive bins is due to the effect that two 
fibres cannot be assigned closer than 55”. 

The two models considered in this section are based on the 
same simulation, assume the same (Planck) cosmology, and have 
identical stellar mass functions. They, nevertheless, predict signifi¬ 
cantly different lensing signals for the most massive LEGs because 
of differing amounts of scatter in their stellar mass - halo mass re¬ 
lations. The main motivation for the present analysis of clustering 
was the hope that it would provide an independent observational 
indication of which model is more realistic. To the extent that it 
provides a significantly better fit to SDSS measurements for both 
kinds of data in both the one-halo and the two-halo regime of the 
two most massive bins, Gil — P' appears to be a higher fidelity 
representation of the real stellar mass - halo mass relation (includ¬ 
ing scatter) than H15f — P'. 


6 LENSING CALIBRATION OF SZ AND X-RAY 
SCALING RELATIONS 

PI3 and A15 stacked the Planck and ROSAT sky maps around 
the LEG sample analysed in this paper in order to measure mean 
Sunyaev-Zeldovich signals and mean X-ray luminosities as func¬ 
tions of stellar mass. The large sample size resulted in detections to 
substantially lower masses and with higher precision than had been 
possible in earlier work. In both papers, G13 — W7' was combined 
with simple assumptions about halo gas profiles and power-law re¬ 
lations between halo mass and observable signal to carry out de¬ 
tailed forward modelling of the stacking and detection procedures. 
This allowed estimation both of effective halo mass as a function of 
LEG stellar mass, and of the parameters of the scaling relation be¬ 
tween observable and halo mass. Since G13 — W7' was originally 
tuned to reproduce the stellar mass function of the SDSS, this pro¬ 
cedure was effectively an abundance-matching calibration of the 
halo mass-stellar mass relation which, as we have seen in earlier 
sections, is model-dependent because of: (i) residual mismatch be¬ 
tween model and observed stellar mass functions; (ii) uncertainties 
in the cosmological parameters; and (iii) the model-dependence of 
the detailed shape of the halo mass distribution at fixed LEG stellar 


mass (including effects due to varying satellite fractions and off¬ 
sets). 

The goal of the current section is to remove most of this 
model-dependence by recalibrating the scaling relations of PI 3 and 
A15 using our gravitational lensing measurements. The new cali¬ 
brations will be model-dependent only to the extent that the ratio 
of the effective halo mass for gravitational lensing to that for SZ 
(or X-ray) stacking varies between models. Since all three effective 
masses correspond to rather similar moments of the halo mass dis¬ 
tribution, this dependence is quite weak, and we include an estimate 
of the associated uncertainty in the analysis leading to our final re¬ 
calibrated scaling relations. When estimating halo masses, P13 and 
A15 matched the observed LEGs to simulated LEGs of the same 
stellar mass in the z = 0 output of G13 — W7'. For consistency, 
we will match in the same way in this section, rather than matching 
in stellar mass, luminosity and redshift, as in earlier sections. 

In the left panel of Fig.[T0]we compare our SDSS measure¬ 
ments of differential lensing surface density profiles with the pre¬ 
dictions of G13 — W7'. Dashed curves are for the original model, 
as used by P13 and A15. These are close to the observations, sug¬ 
gesting that the model already represents the halos of SDSS galax¬ 
ies well, and that the shifts needed to recalibrate to the lensing 
data will be small. For each stellar mass bin, we estimate the shift 
(and its uncertainty) by assuming that the halos of all model LEGs 
change mass by A log M 200 without changing concentratiorQ. In 
the one-halo regime, the differential density profiles then shift up¬ 
wards by IA log M 200 and to the right by the same amount. We 
find the best shift and its uncertainty by calculating values for 
the model curves with respect to the SDSS data at Vp < 4i7200, the 
region over which iHavashi & Whitel bOOSh found the lensing pro¬ 
files of central galaxies to be well described by the one-halo term. 
For the three most massive stellar mass bins, we exclude the inner¬ 
most three, two or one SDSS data points from the fit (indicated by 
a lighter colour using upper triangle symbols in Fig. [Jo) because 
measurements on these scales may be affected by light from the 
central galaxy. We determi ne appropriate inner scales for exclusion 
using the criteria set out bv iMandelbaum et alj ( l2005h . 

The solid curves in the left panel of Fig. [TO] show results for 
G13 — W7' after applying these recalibration shifts. A short verti¬ 
cal bar on each curve marks Vp = 4i?200- The shifts themselves are 
shown as the black points in the right panel, with the shorter error 
bars indicating uncertainties from the fitting. Note the quantity 
plotted for the x-axis, log(M*), is averaged by taking equation |4] 
for each galaxy as weights. The SDSS measurements for the lowest 
stellar mass bin are very noisy, but this point is not needed for our 
recalibration below. At higher mass, both the shifts and their statis¬ 
tical uncertainties are quite small. Assuming that the shape of the 
halo mass distribution at given stellar mass is unchangec[3, these 
shifts can be applied to the effective halo masses estimated by P13 
and A15 in order to calibrate them to the lensing data. To obtain 
the overall uncertainties shown in Fig. [TO] we add in quadrature a 
systematic uncertainty due to the fact that the ratio of lensing to 
SZ (or X-ray) effective mass may differ between G13 — W7' and 
the real galaxies. This is estimated as the larger of columns 3 and 

^ A log M 200 is fairly small for G 13 — W7'. Given the flat slope of con- 
centration halo mass relations te.g. lBullock et aflboOlUDuffy et al.ll2008l 
where the slope is about -0.1). the assumption of unchanged concentration 
should be reasonable, i.e., 0.1 dex of change in halo masses corresponds to 
at most a few percent change in halo concentrations. 

® This should be a reasonable assumption given the fact that the curve in 
the right panel is close to flat at log(M*) > 10.5 
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Figure 10. Left: As Fig.|4] but now comparing our SDSS lensing measurements {dots and triangles with errors) to the predictions of G13 — W7', the model 
used by P13 and A15. Dashed cuiwes give predictions for the original model, while solid curves are shifted upwards and to the right (or down and to the left) 
by (A log M2 oo)/ 3, the (potentially negative) amount needed to bring them into optimal agreement with the SDSS data in the one-halo regime, rp < 4i7200 
(indicated by a short vertical bar on each curve). Upper triangles are shown in a lighter colour and are not used when fitting because of concerns that they 
may be affected by light from the LEG. To be consistent with P13 and A15, each LEG in the SDSS was matched in stellar mass alone to LEGs in the 
z = 0 output of the simulation when making this plot. Right: Recalibration shifts in log M 200 ^fe shown as filled symbols with en'or bai's indicating: (a) 
statistical uncertainties in the shifts resulting from the error bai's on the SDSS lensing measurements; and (b) systematic uncertainties (taken from TableO in 
the conversion from effective lensing mass to (log M 200 ). The shorter eri'or bai' reflects (a) alone, whereas the longer includes (a) and (b) added in quadrature. 
The pink shaded region interpolates linearly both in the recalibration shift and in its uncertainty. 


5 in Tabled which give the scatter in similar effective mass ratios 
over our 16 galaxy formation models. The black line and the sur¬ 
rounding red band in Fig.^Olthen give the correction and associated 
uncertainty which we apply to the effective halo masses of P13 and 
A15 in order to recalibrate their scaling relations. Note that the un¬ 
certainties come predominantly from noise in the lensing results at 
small stellar mass, but from modelling systematics at large stellar 
mass. 

In Fig. we show the effect of recalibration on the scaling 
relations. The upper left panel replots the Vsoo-ATsoo relation from 
Fig. 9 of P13 as black symbols. As in the original paper, these have 
vertical error bars, indicating flux uncertainties estimated from the 
stacking analysis, but allow for no uncertainty in effective halo 
mass Msoo- Red symbols show the recalibrated relation. The verti¬ 
cal error bar is unchanged but the points now also have a horizontal 
error baS corresponding to the pink band in Fig. [10] The upper 
right panel is a similar plot for the I/totai-ATsoo relation from Ta¬ 
bles 3 and A1 of A15. Exact definitions of the various observational 
quantities plotted here can be found in the original papers. We fit 
these relations to power-law expressions similar to those of Eqn. 1 
in PI3 and Eqn. 3 in A15. For the original data we minimise a 


® Uncertainties for the twelve SZ/X-ray stellar mass bins are estimated 
by interpolating linearly between those for the neighbouring lensing bins. 
Since we have changed to the Planck cosmology, which increases the stellar 
mass of each SDSS galaxy by about 10% over that used by P13 and A15, an 
additional con'ection of log(0.704/0.673)^ was added to the central stellar 
mass of each SZ/X-ray bins bin before inteipolating. 


estimated from the vertical offsets in ysoo and Ltotai and the boot¬ 
strap errors from the original papers. For the recalibrated data, we 
take account of the uncertainties in both directions. This is a com¬ 
plex statistical problem because of the strongly correlated errors 
on neighbouring bins resulting from the fact that our twelve cor¬ 
rection factors are interpolated/extrapolated from the seven lensing 
measurements of Fig. El This effect is particularly important for 
the highest stellar mass bin s. We have d eveloped a fully Bayesian 
procedure similar to that of lKellvl ( l2007h and presented in detail in 
the Appendix, which accounts consistently for these effects. 

We write our scaling relations as 

and 

Ltotai = To X E{zy^^ X (8) 

These differ from the scaling relations of PI3 and A15 only in 
that we here adopt lower mass pivot-points than in the original pa¬ 
pers in order to decorrelate uncertainties on the two parameters. 
This is necessary because our lensing data constrain halo masses 
most tightly at a few times 10^® M©. In Eqn. [8] is a bolo- 

metric correction factor to convert the 0.5- 2.0 keV luminosities t o 
total luminosities, and is given in Table 2 of lAnderson et al.l ( l2015h . 
E(z) = -f FIa is the dimensionless Hubble param¬ 

eter in a flat cosmological model to allow for self-similar evolution, 
which is already included in the P13 definition of Y 50 Q (see the end 
of their Section 1). 
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Figure 11. The SZ I 500 -AT 500 relation (top left) and X-ray Ltotal-A^soo relation (top right). Black dots with vertical en'or bars are the original relations 
given in P13 and A15. Red dots with both vertical and horizontal error bars are the same measurements after recalibrating the effective halo masses using our 
lensing measurements. The black and red dashed lines are best-fit power-law scaling relations for the original and the recalibrated data, respectively. Triangles 
along the red dashed lines indicate the maximum likelihood effective halo masses associated with each stellar mass bin. The bolometric correction, is 

included for the X-ray model. This varies between bins, so the model is not strictly a straight line. At each halo mass, the pink shaded regions enclose 68% 
of the scaling relations proposed by our MCMC chains. The bottom panels show the best-fit parameters from these chains as symbols surrounded by contours 
indicating their 1- and 2-a confidence regions. Black and red again indicate the original and recalibrated cases, respectively. The confidence regions are larger 
after recalibration because both systematic modelling uncertainties and random noise in the effective halo mass estimates are now consistently included. 


The best-fit scaling relations are overplotted in the upper pan¬ 
els of Fig.[TT]as black dashed lines for the original measurements, 
and red dashed lines for the recalibrated measurements. The lower 
panels show the amplitudes and slopes of these fits, together with 
contours indicating 1 and 2-0 confidence regions. For the origi¬ 
nal data, these are ellipses obtained from the derivatives of at its 
minimum. For the recalibrated data, they are taken directly from the 
MCMC chains used to map out the likelihood surfaces. Red trian¬ 
gles along the best fit relations show the maximum likelihood effec¬ 
tive halo masses associated with each stellar mass bin. At each halo 
mass, the pink shaded regions enclose 68% of the scaling relations 
proposed by our MCMC chains. The maximum likelihood parame¬ 
ters for our recalibrated scaling relations and their 68% confidence 
regions, both obtained from the one-dimensional marginalised like¬ 
lihood distributions generated by our MCMC chains, are: 

10® X Fm = 2.31 [1.93, 2.67], om = 1.61 [1.43,1.75] (9) 
and 

Lo/10'‘^erg s"^ = 1.34 [1.12,1.67], a = 1.91 [1.78, 2.07] 

( 10 ) 

The recalibration makes the slope of the best-fit SZ relation 


slightly shallower, while its amplitude increases by about 30%. In 
addition, the confidence region is significantly larger now that un¬ 
certainties in the estimation of effective halo masses are consis¬ 
tently included. For the best fit X-ray scaling relation, the slope 
barely changes, but the amplitude increases by about 40%. These 
changes in amplitude are due primarily to the decrease in effective 
halo mass for most stellar mass bins (corresponding to the negative 
corrections in Fig.llOt. Note that, as expected, the uncertainties in 
the amplitude and slope of the scaling relations are only weakly 
correlated, reflecting the fact that mean SZ and X-ray fluxes are 
best constrained (to better than ±20% at one o) for halo masses 
near our new pivot value. This is perhaps the most interesting result 
of our paper, since the systematic X-ray and SZ properties of halos 
in the range 10 ^^ ® < M 500 /Mq < lO^'^ have not previously been 
constrained at anything approaching this level of precision because 
of difficulties in the selection of representative samples of objects, 
in the estimation of their halo masses, and in the measurement of 
their X-ray and SZ signals. 
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7 CONCLUSIONS 


We have measured weak gravitational lensing profiles at high 
signal-to-noise and as a function of stellar mass for stacks of Lo¬ 
cally Brightest Galaxies selected from SDSS/DR7, comparing the 
observed signal over the projected radius range 30h“^kpc < Vp < 
30h”^Mpc to predictions from eight semi-analytic galaxy forma¬ 
tion simulations differing in underlying N-body simulation, in as¬ 
sumed cosmological parameters and in the modelling of galaxy 
formation processes. This L BG sample was previously used by 
PlanckCollalroratioj^( |20 1 3h to measure stacked SZ signal and by 
Anderson et al. l l2015h to measure stacked X-ray signal as functions 
of LBG stellar mass. Both studies used a galaxy formation simula¬ 
tion from the set in this paper to carry out forward modelling of the 
sample selection and signal measurement processes, and hence to 
derive scaling relations between the mass of dark matter halos and 
the properties of their hot gas atmospheres. Here, we use our sim¬ 
ulation set to explore the model-dependence of such calibrations, 
before recalibrating the SZ and X-ray scaling relations to a halo 
mass scale based directly on our lensing results. This recalibration 
takes full account of the observational uncertainties in the lensing, 
SZ and X-ray observations, as well as residual modelling uncer¬ 
tainties arising from the scatter in halo mass at given LBG stellar 
mass. 

We compared the differential surface density profiles mea¬ 
sured around our stacks of SDSS LBGs to those around LBGs de¬ 
fined in an exactly analogous way in our galaxy formation simula¬ 
tions, taking care to reproduce both the redshift distribution of the 
observed galaxies at each stellar mass and the redshift-dependent 
weighting of their lensing signals. For each of our eight simula¬ 
tions we compared to predictions from the simulation both in its 
original form and after applying small stellar mass corrections to 
bring its stellar mass function into exact agreement with that ob¬ 
served in the SDSS. Such corrections eliminate inaccurate stellar 
mass functions as a possible source of disagreement between the 
simulated and observed lensing signals. From these comparisons 
we found: 


• All of our galaxy formation simulations predict lensing pro¬ 
files that are in good qualitative agreement with the SDSS data 
over two orders of magnitude in projected radius and a factor 
of 30 in stellar mass. The quantitative level of agreement varies 
between simulations, however, with amplitude offsets of up to a 
few tenths of a dex i n som e cas es. The specific model used by 
IPlanck CollaborationI ( 1201 3h and lAnderson et al.l l l2015h to cali¬ 
brate their scaling relations produces one of the best fits to the lens¬ 
ing data and requires very small corrections to bring its stellar mass 
function into agreement with SDSS. 

• Including stellar mass corrections improves the agreement in 
mean halo mass as a function of stellar mass between simulations, 
but significant scatter remains, particularly for large stellar masses. 
In most cases, these corrections also improve agreement with the 
SDSS lensing measurements for logM, < 10 ^^’^Mq, but this is 
often not the case at high stellar mass. 

• Models with identical stellar mass functions (as a result of the 
stellar mass corrections) and based on the same N-body simulation 
can produce significantly different lensing profiles and halo mass 
distributions at fixed stellar mass. This is especially noticeable for 
logM* > 1O^^'"^M0. In general, shifts in the amplitude of the 
lensing profiles are matched, at least qualitatively, by shifts in mean 
host halo mass. Because the scatter in halo mass is large and model- 
dependent, LBG stellar mass has large systematic uncertainties as 


a proxy for host halo mass, and so may not be useful for many 
cosmological applications. 

• Simulations with the same stellar mass function, the same 
galaxy formation model, but different cosmologies can make dif¬ 
ferent predictions both for lensing profiles and for halo mass dis¬ 
tributions. The differences are not straightforward to interpret, with 
several factors playing a significant role. Again, this is likely to 
complicate the drawing of conclusions about cosmological param¬ 
eters from lensing observations of the type considered here. 

Because the model-dependence of the scatter in halo mass at 
fixed LBG stellar mass leads to systematic uncertainties in the in¬ 
terpretation of the observed SDSS lensing profiles, we investigated 
whether additional useful constraints could be obtained from the 
clustering of other galaxies around LBGs, specifically a reference 
sample of all objects more massive than logM* > 10^°M©. A 
comparison of predictions from two of our simulations to cluster¬ 
ing measurements from SDSS showed that differences in effective 
LBG halo mass inferred from lensing are also seen as differences 
in the amplitude of the large-scale (two-halo) clustering of other 
galaxies around LBGs, and that conclusions about the relative mer¬ 
its of the models drawn from their lensing predictions are consis¬ 
tent with those drawn from their large-scale clustering. On smaller 
scales, clustering is sensitive to how reference galaxies populate 
LBG halos, and this introduces additional model-dependencies. In 
our case, models which predict well the observed mean distribution 
of dark matter around LBGs on scales of a few hundred kpc do not, 
for LBGs of intermediate stellar mass, predict well the distribution 
of reference galaxies on these same scales. 

By combining results from our eight simulations both with 
and without stellar mass corrections we are able to estimate, for 
each LBG stellar mass bin in our lensing analysis, the system¬ 
atic uncertainty in the ratio between the effective halo mass for 
gravitational lensing and the mean and median halo masses. Such 
uncertainties reflect the model-dependence of the detailed shape 
of the halo mass distribution at given LBG stellar mass. In Sec¬ 
tion 6 we use them to estimate the systematic uncertainty in the 
ratio of effective halo mass for gravitational lensing to effective 
halo mass for S Z or X-ray stacking when r ecali brating the scal¬ 
ing re lations of IPlanck CollaborationI (1201 3h and lAnderson et all 
( |2015H^°L Because the lensing analysis is carried out in broader stel¬ 
lar mass bins than the earlier stacking analyses, interpolation of the 
lensing results is required when carrying through the recalibration. 
This complicates the formal problem of accounting properly both 
for observational uncertainties in the three stacking measurements 
and for systematic uncertainties from residual model dependencies. 
We have developed a Bayesian analysis technique for this problem 
which is presented in the Appendix. We believe the resulting recali¬ 
brated scaling relations, given in equations 8 to 10, to be robust and 
minimally model dependent and to be quoted with realistic uncer¬ 
tainties. They should be fully representative for halos in the mass 
range ® < Msoo/L'^Mq < 10^^'®, a significantly broader 
and lower range than it has been possible to cover through obser¬ 
vations of individual objects. 


P13 and A15 have made other assumptions which may introduce sys¬ 
tematic uncertainties. For the SZ scaling relation, a pressure profile has to 
be assumed to estimate the effective halo mass. To obtain the X-ray scal¬ 
ing relation, A15 have made bolometric corrections to the observed X-ray 
flux, which might introduce additional uncertainties. Studying such effects 
is outside the scope of this paper. 
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APPENDIX 

This appendix describes how we fit the power-law scaling relations 
of Equations |7] and [8] to the recalibrated data shown as red points 
with both horizontal and vertical error bars in the upper panels of 
EigdU This fitting problem is difficult because of the need to deal 
with uncertainties on both variables which vary from point to point, 
and are strongly correlated between neighbouring points. Since the 
recalibrated scaling relations are the principal quantitative result 
of our paper, we considered it important to devise a well defined 
and rigorous statistical procedure to derive them. Eor simplicity, 
the following exposition considers the SZ case and refers to the 
two variables as Y and M rather than as Ysoo and M 500 . The X- 
ray case can be treated in an exactly ana logous way. T he method 
we propose is close to those developed bv iKellvl ll2007h . 

The data to be fitted consist of twelve estimates, (Yi, Mi), one 
for each of the twelve stellar mass bins we consider. The Yi values 
come directly from PI3, whereas the P13 Mi values are recali¬ 
brated using the correction indicated by a solid black line in Pig. 


(right hand panel). Let us denote as yi and rrn the true mean SZ 
flux and the true effective halo mass for the population of galaxies 
with stellar masses in the range of bin i. Our stochastic model may 
be written as 

Yi = yi + Si (A-1) 

and 

log Mi = log rrii + ei + r]i, (A-2) 

where Si, ti and ry are zero-mean, Gaussian deviates. The first. 
Si, reflects the observational uncertainty in the stacked SZ flux that 
PI3 measured for bin i, and should be uncorrelated both between 
bins and with uncertainties in effective halo mass. Hence 

(SiSj) = Cy,ij, (Sitj) = (Sirjj) = 0, (A-3) 

where the covariance matrix Cy,ij is diagonal with elements given 
by the (square of the) bootstrap errors estimated by PI3. We take 
the distribution of SZ flux to be Gaussian in Y rather than in log Y 
in order to deal with lower signal-to-noise bins where the estimated 
flux turns out to be negative in one case. 

In contrast, for the effective mass it is more natural to assume 
the uncertainties to be Gaussian in log M. We distinguish two noise 
sources. The larger, Si, reflects uncertainty in the correction from 
the effective mass estimated by PI3 to that implied by our lensing 
measurements. The twelve ti are thus linearly related to the seven 
6 a which represent the differences between the seven correction 
estimates of Pig.[T0]and their true population values. These 6 a rep¬ 
resent both systematic modelling uncertainties and observational 
noise in the lensing data. They can be taken to be independent, 
zero-mean, Gaussian variates with variances equal to the (squares 
of the) error bars in Pig.[T0l However, the ti = Diada will have a 
nondiagonal (and indeed, singular) covariance matrix. 

Pinally, the rji are introduced as a mathematical device to 
make Cm,ij well behaved and invertible. Additionally, the r]i can 
be thought of as accounting for errors resulting from the interpo¬ 
lation of correction factors derived for our relatively broad lensing 
stellar mass bins to the narrower stellar mass bins used by P13 and 
A15. This is expected to induce relatively small additional errors, 
so for simplicity, we just take all the yi to be independent, zero- 
mean Gaussian variates with rms a = 0.03dex. With these assump¬ 
tions we have 

(eiEj) = Cm.ij, {rjiyj) = a'^hj, {€ir]j)=0, (A-4) 

where lij is the identity matrix. We can then define a nonsingular 
covariance matrix 

= ((log Mi - log rrii) (log Mj-log rrij)) = Cm.ij+a^Iij- 

(A-5) 

Eor a power-law scaling relation of the form 

log y — a log m -I- ti, (A-6) 

we can now write the likelihood of a specific model (specified by 
a, b and the twelve rrii) given the data [the twelve (Mi, Y)] as that 
corresponding to a multivariate Gaussian 

2 log £ = -(log Mi - log mi)C'~}j (log Mj - log nij) (A-7) 
-{Yi - yi)C-UY^ - Vi), 

where log yi = a log mi — b, the rrii can be thought of as nuisance 
parameters, and we have dropped a constant which is independent 
of the model parameters. Note that it was important to include the 
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additional stochastic element represented by the r)i in order to en¬ 
sure that the inverse covariance matrix in the first term is well de¬ 
fined. This likelihood function can be fed to a Monte Carlo Markov 
Chain program in order to explore the likelihood surface, and, in 
particular, to project it onto the (a, b) plane where we can iden¬ 
tify a best-fit (i.e. maximum likelihood) scaling relation and define 
confidence intervals for its parameters. 
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